In [62]:

    

%pylab inline
from scipy import signal
import math
pylab.rcParams['figure.figsize'] = (16, 6)
pylab.rcParams["font.size"] = "20"


    

    




Populating the interactive namespace from numpy and matplotlib

In [91]:

    

amp = 0.5      # amplitude of the cosine wave
fc1 = 8.0      # frequency of the cosine wave
fc2 = 48.0
phase = 30.0   # desired phase shift of the cosine in degrees
fs = 100.0     # sampling frequency, 100 samples/s
ts = 1.0 / fs  # sample interval
N = 256        # FFT size
tmax = (2.0*N)/fs # length of signal in seconds
t = arange(0.0, tmax, ts)


    

In [92]:

    

phi = deg2rad(phase) # phase shift to radians
# time domain signal with phase shift
sig = amp * cos(2.0 * pi * fc1 * t + phi) + \
      amp/4 * cos(2.0 * pi * fc2 * t + phi)
# and plot
grid(True)
plot(t,sig)
title('%3.1f cos(2 pi %d t + pi/6)' % (amp, fc1))
xlabel('time (sec)')
ylabel('x(t)')


    

    
Out[92]:





<matplotlib.text.Text at 0x7f7d6c0e5650>






    

In [93]:

    

X = fft.rfft(sig) / float(N) # N-point real DFT
frq = fft.rfftfreq(sig.size, d=ts) # get frequency bins


    

In [94]:

    

# verify things
print X.size, sig.size, t.size, frq.size, N
print "t  :", t[0], t[-1] 
print "frq:", frq[0], frq[-1]
peak_ix = argmax(abs(X))
print "peak of %f at ix %d, %f Hz" % (abs(X[peak_ix]), peak_ix, frq[peak_ix])


    

    




257 512 512 257 256
t  : 0.0 5.11
frq: 0.0 50.0
peak of 0.498500 at ix 41, 8.007812 Hz

In [95]:

    

grid(True)
stem(frq, abs(X)) # magnitudes vs frequencies
xlabel('f (Hz)')
ylabel('|X(k)|')


    

    
Out[95]:





<matplotlib.text.Text at 0x7f7d6c209bd0>






    

In [97]:

    

mx = max(abs(X))
threshold = mx/10  # tolerance threshold
# maskout values that are below the threshold
X2 = [x if abs(x)>threshold else 0+0j for x in X]
phase = rad2deg(arctan2(imag(X2),real(X2))) # phase information
grid(True)
stem(frq,phase) # phase vs frequencies
xlabel('f (Hz)')
ylabel('dX(k)')


    

    
Out[97]:





<matplotlib.text.Text at 0x7f7d6b2b85d0>






    

In [ ]:

    

sig_recon = N * fft.irfft(X) # reconstructed signal
grid(True)
plot(t,sig)
title('reconstructed')
xlabel('time (sec)')
ylabel('x(t)')


    

In [ ]: